#################################################################################
### Racialized Misinformation, Factual Corrections, and Prejudicial Attitudes ###
### [Main Replication Code for Study 3]                                       ###
### Authors: Eddy S. F. Yeung, Joseph Glasgow                                 ###
### Date: June 26, 2025                                                       ###
#################################################################################

### Set-up ----
## Clean the working environment and set the working directory
rm(list = ls())
setwd("~/Desktop/racialized_misinfo/replication/Study 3") # set your working directory here, which should also contain the survey data ("study3_dataset.csv")

## Import the dataset
df <- read.csv("study3_dataset.csv")

## Load required packages
library(tidyverse)
library(texreg)
library(estimatr)
library(cowplot)
library(grid)
library(gridExtra)
library(interflex)
library(multcomp)
library(broom)
library(ivreg)
library(modelsummary)
library(parameters)
library(effectsize)

### Code demographic variables ----
## Age
df$age <- df$yob + 12

## Female (= 1)
df$female <- df$gender - 1

## Race (0 = White; 1 = Black; 2 = Hispanic; 3 = other)
df <- df %>% mutate(race4 = case_when(
  racial == 1 ~ 0,
  racial == 2 ~ 1,
  racial == 3 ~ 2,
  racial >= 4 ~ 3
))

## White (= 1)
df$white <- ifelse(df$racial == 1, 1, 0)

## Married (= 1)
df$married <- ifelse(df$marital == 1, 1, 0)

## Education (less than high school = 0; doctorate degree = 6)
df <- df %>% mutate(education = case_when(
  edu == 1 ~ 0,
  edu == 2 ~ 1,
  edu == 3 ~ 2,
  edu == 4 ~ 3,
  edu == 5 ~ 4,
  edu == 6 | edu == 7 | edu == 8 ~ 5,
))

## Income (less than $10,000 = 0; $150,000 or more = 12)
df <- df %>% mutate(income = case_when(
  income == 1 ~ 0,
  income == 2 | income == 3 ~ 1,
  income == 4 | income == 5 ~ 2,
  income == 6 ~ 3,
  income == 7 ~ 4,
  income == 8 ~ 5,
  income == 9 ~ 6,
  income == 10 ~ 7,
  income == 11 ~ 8,
  income == 12 ~ 9,
  income == 13 ~ 10,
  income == 14 ~ 11,
  income == 15 | income == 16 | income == 17 ~ 12
))

## Party identification (0 = strong Democrat; 6 = strong Republican)
df <- df %>% mutate(pid = case_when(
  pid1 == 1 & pid_d == 1 ~ 0,
  pid1 == 1 & pid_d == 2 ~ 1,
  (pid1 == 3 | pid1 == 4) & pid_i == 2 ~ 2,
  (pid1 == 3 | pid1 == 4) & pid_i == 3 ~ 3,
  (pid1 == 3 | pid1 == 4) & pid_i == 1 ~ 4,
  pid1 == 2 & pid_r == 2 ~ 5,
  pid1 == 2 & pid_r == 1 ~ 6
))

## Democrat (= 1)
df$dem <- ifelse(df$pid == 0 | df$pid == 1 | df$pid == 2, 1, 0)

## Republican (= 1)
df$gop <- ifelse(df$pid == 4 | df$pid == 5 | df$pid == 6, 1, 0)

## Self-reported ideology (0 = extremely liberal; 6 = extremely conservative)
df$ideo <- df$leftright - 1

## Attentiveness (0 = lowest; 2 = highest)
df$screener_2[is.na(df$screener_2)] <- 0
df$screener_2_correct <- ifelse(df$screener_2 == 3, 1, 0)
df$screener_4[is.na(df$screener_4)] <- 0
df$screener_4_correct <- ifelse(df$screener_4 == 2, 1, 0)
df$attentive <- df$screener_2_correct + df$screener_4_correct

## Impute preregistered covariates
df$age[is.na(df$age)] <- mean(df$age, na.rm = TRUE)
df$female[is.na(df$female)] <- mean(df$female, na.rm = TRUE)
df$white[is.na(df$white)] <- mean(df$white, na.rm = TRUE)
df$married[is.na(df$married)] <- mean(df$married, na.rm = TRUE)
df$education[is.na(df$education)] <- mean(df$education, na.rm = TRUE)
df$income[is.na(df$income)] <- mean(df$income, na.rm = TRUE)
df$ideo[is.na(df$ideo)] <- mean(df$ideo, na.rm = TRUE)
df$pid[is.na(df$pid)] <- mean(df$pid, na.rm = TRUE)

### Code outcome variables ----
## Affective prejudice (0 = most favorable and warm; 100 = least favorable and warm)
df$affective_prej <- 100 - df$affect_1

## Desire for social distance (0 = min; 100 = max)
df$social_distance <- 15 - (df$social_distance_1 + df$social_distance_2 + df$social_distance_3)
df$behavioral_prej <- df$social_distance / 3 * 25

## Support for immigration (0 = decrease a lot; 4 = increase a lot)
df$immi_support <- 5 - df$immi_support

## Direct manipulation check (1 = correct)
df$mani_check <- df$fmc - 1

## Perceived effect of illegal immigration on crime rate (0 = decrease a lot; 6 = increase a lot)
df <- df %>% mutate(perceived_crime = case_when(
  perceived_crime_2b == 1 ~ 0,
  perceived_crime_2b == 2 ~ 1,
  perceived_crime_2b == 3 ~ 2,
  perceived_crime_1  == 3 ~ 3,
  perceived_crime_2a == 3 ~ 4,
  perceived_crime_2a == 2 ~ 5,
  perceived_crime_2a == 1 ~ 6
))

## Perceived credibility of the vignette (0 = very untrustworthy; 4 = very trustworthy)
df$cred <- 5 - df$cred

## Overestimation of the crime rate of undocumented immigrants (truth = 4.2)
df$overest_UI <- df$pre_est_3 - 4.2

## Overestimation of the crime rate of legal immigrants (truth = 8.1)
df$overest_LI <- df$pre_est_2 - 8.1

### Figure 5 ----
## Average level of misperception about undocumented immigrants in the full sample (true = 4.2)
mean(df$pre_est_3, na.rm = T)
median(df$pre_est_3, na.rm = T)
sum(!is.na(df$pre_est_3))

## Pretreatment estimate of the crime rate of undocumented immigrants by party ID
group_mean_temp1 <- df %>% 
  filter(gop == 0 | gop == 1) %>% 
  group_by(gop) %>% 
  summarize(xvalue = mean(pre_est_3, na.rm = T))
p0.1 <- 
  ggplot(data = subset(df, gop == 0 | gop == 1), 
         aes(x = pre_est_3, group = as.factor(gop), fill = as.factor(gop))) +
  geom_density(alpha = .25) +
  geom_vline(data = group_mean_temp1,
             aes(xintercept = xvalue, color = as.factor(gop)), 
             linetype = "dashed", linewidth = 0.5) +
  scale_color_manual(labels = c("Non-GOP", "GOP"), values = c("#00008B", "#EFC000FF")) +
  scale_fill_manual(labels = c("Non-GOP", "GOP"), values = c("#00008B", "#EFC000FF")) +
  geom_vline(xintercept = 4.2, color = "red", linetype = "solid", linewidth = 1) +
  xlab("Estimated Crime Rate of Undocumented Immigrants") +
  ylab("Density") +
  ggtitle("Prior Beliefs about Undocumented Immigrants\nby Partisanship") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 11),
        axis.text = element_text(color = "black", size = 11),
        legend.justification = c(1, 1), legend.position = c(.95, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(.75, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.45, size = 13, face = "italic"))

## Pretreatment estimate of the crime rate of undocumented immigrants by racial identity
group_mean_temp2 <- df %>% 
  filter(white == 0 | white == 1) %>% 
  group_by(white) %>% 
  summarize(xvalue = mean(pre_est_3, na.rm = T))
p0.2 <- 
  ggplot(data = subset(df, white == 0 | white == 1), 
         aes(x = pre_est_3, group = as.factor(white), fill = as.factor(white))) +
  geom_density(alpha = .25) +
  geom_vline(data = group_mean_temp2,
             aes(xintercept = xvalue, color = as.factor(white)), 
             linetype = "dashed", linewidth = 0.5) +
  scale_color_manual(labels = c("Nonwhite", "White"), values = c("#00008B", "#EFC000FF")) +
  scale_fill_manual(labels = c("Nonwhite", "White"), values = c("#00008B", "#EFC000FF")) +
  geom_vline(xintercept = 4.2, color = "red", linetype = "solid", linewidth = 1) +
  xlab("Estimated Crime Rate of Undocumented Immigrants") +
  ylab("Density") +
  ggtitle("Prior Beliefs about Undocumented Immigrants\nby Racial Identity") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 11),
        axis.text = element_text(color = "black", size = 11),
        legend.justification = c(1, 1), legend.position = c(.95, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(.75, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.45, size = 13, face = "italic"))

## Average level of misperception about legal immigrants in the full sample (true = 8.1)
mean(df$pre_est_2, na.rm = T)
median(df$pre_est_2, na.rm = T)
sum(!is.na(df$pre_est_2))

## Pretreatment estimate of the crime rate of undocumented immigrants by party ID
group_mean_temp3 <- df %>% 
  filter(gop == 0 | gop == 1) %>% 
  group_by(gop) %>% 
  summarize(xvalue = mean(pre_est_2, na.rm = T))
p0.3 <- 
  ggplot(data = subset(df, gop == 0 | gop == 1), 
         aes(x = pre_est_2, group = as.factor(gop), fill = as.factor(gop))) +
  geom_density(alpha = .25) +
  geom_vline(data = group_mean_temp3,
             aes(xintercept = xvalue, color = as.factor(gop)), 
             linetype = "dashed", linewidth = 0.5) +
  scale_color_manual(labels = c("Non-GOP", "GOP"), values = c("#00008B", "#EFC000FF")) +
  scale_fill_manual(labels = c("Non-GOP", "GOP"), values = c("#00008B", "#EFC000FF")) +
  geom_vline(xintercept = 8.1, color = "red", linetype = "solid", linewidth = 1) +
  xlab("Estimated Crime Rate of Legal Immigrants") +
  ylab("Density") +
  ggtitle("Prior Beliefs about Legal Immigrants\nby Partisanship") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 11),
        axis.text = element_text(color = "black", size = 11),
        legend.justification = c(1, 1), legend.position = c(.95, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(.75, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic"))

## Pretreatment estimate of the crime rate of undocumented immigrants by racial identity
group_mean_temp4 <- df %>% 
  filter(white == 0 | white == 1) %>% 
  group_by(white) %>% 
  summarize(xvalue = mean(pre_est_2, na.rm = T))
p0.4 <- 
  ggplot(data = subset(df, white == 0 | white == 1), 
         aes(x = pre_est_2, group = as.factor(white), fill = as.factor(white))) +
  geom_density(alpha = .25) +
  geom_vline(data = group_mean_temp4,
             aes(xintercept = xvalue, color = as.factor(white)), 
             linetype = "dashed", linewidth = 0.5) +
  scale_color_manual(labels = c("Nonwhite", "White"), values = c("#00008B", "#EFC000FF")) +
  scale_fill_manual(labels = c("Nonwhite", "White"), values = c("#00008B", "#EFC000FF")) +
  geom_vline(xintercept = 8.1, color = "red", linetype = "solid", linewidth = 1) +
  xlab("Estimated Crime Rate of Legal Immigrants") +
  ylab("Density") +
  ggtitle("Prior Beliefs about Legal Immigrants\nby Racial Identity") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 11),
        axis.text = element_text(color = "black", size = 11),
        legend.justification = c(1, 1), legend.position = c(.95, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(.75, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic"))
prior <- plot_grid(p0.3, p0.4, p0.1, p0.2, labels = "AUTO", nrow = 2, label_fontfamily = "Times")
ggsave(file = "Figure 5.pdf", prior, width = 8, height = 6)
# ggsave(file = "Figure 5.tif", prior, dpi = 300, width = 8, height = 6)

### Figure 6 and Figure S16----
df$treatment <- as.factor(df$immi_treatment)

## Affective prejudice
# Estimate the ATE
aff_prej_ATE <- lm_robust(affective_prej ~ treatment + age + female + white + 
                            married + education + income + ideo + pid,
                          data = df)

# Estimate the effect size of the ATE
effectsize::cohens_d(affective_prej ~ treatment, data = df)

# Visualize the results
group_mean_temp1 <- df %>% 
  group_by(treatment) %>% 
  summarize(xvalue = mean(affective_prej, na.rm = T))
p1 <- 
  ggplot(data = df, aes(x = affective_prej, group = treatment, fill = treatment)) +
  geom_density(alpha = .25) +
  geom_vline(data = group_mean_temp1,
             aes(xintercept = xvalue, color = treatment), 
             linetype = "dashed", linewidth = 0.6) +
  scale_color_manual(labels = c("Control", "Treatment"), values = c("#00008B", "#EFC000FF")) +
  scale_fill_manual(labels = c("Control", "Treatment"), values = c("#00008B", "#EFC000FF")) +
  xlab("Affective Prejudice against Immigrants\n(0 = lowest; 100 = highest)") +
  ylab("Density") +
  ggtitle("Affective Prejudice in Full Sample") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.95, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic"))
p1 <- p1 + 
  annotate(geom = "text", x = 75, y = 0.011, 
           label = paste("beta == ", round(aff_prej_ATE$coefficients[2], 2)),
           size = 4, family = "Times", parse = T) +
  annotate(geom = "text", x = 75, y = 0.010, 
           label = paste("(SE == ", round(aff_prej_ATE$std.error[2], 2), ")"),
           size = 4, family = "Times", parse = T)

# CATE by party ID
lm_robust(affective_prej ~ treatment + age + female + white + married +
            education + income + ideo,
          data = df,
          subset = (dem == 1))
effectsize::cohens_d(affective_prej ~ treatment, data = subset(df, dem == 1))

lm_robust(affective_prej ~ treatment + age + female + white + married +
            education + income + ideo,
          data = df,
          subset = (dem == 0))
effectsize::cohens_d(affective_prej ~ treatment, data = subset(df, dem == 0))

# HTE by seven-point party ID
out <- interflex(Y = "affective_prej", D = "treatment", X = "pid",
                 Z = c("age", "female", "white", "married", "education", "income", "ideo"), 
                 data = df, na.rm = TRUE,
                 estimator = "binning", vcov.type = "robust")
p1_HTE <- 
  plot(out, theme.bw = TRUE, cex.axis = 0.8, cex.lab = 0.8, bin.labs = TRUE,
       xlab = "Party Identification\n(0 = strong Democrat; 6 = strong Republican)", 
       ylab = "Marginal Effect on Affective Prejudice",
       ylim = c(-15, 15)) +
  ggtitle("HTE on Affective Prejudice by Party Identification") + 
  theme(plot.title = element_text(hjust = 0.85, size = 12, face = "italic", family = "Times"))

# CATE by racial identity
lm_robust(affective_prej ~ treatment + age + female + married +
            education + income + ideo + pid,
          data = df,
          subset = (white == 1))
effectsize::cohens_d(affective_prej ~ treatment, data = subset(df, white == 1))

lm_robust(affective_prej ~ treatment + age + female + married +
            education + income + ideo + pid,
          data = df,
          subset = (white == 0))
effectsize::cohens_d(affective_prej ~ treatment, data = subset(df, white == 0))

## Desire for social distance
# Estimate the ATE
beh_prej_ATE <- lm_robust(behavioral_prej ~ treatment + age + female + white + 
                            married + education + income + ideo + pid,
                          data = df)

# Estimate the effect size of the ATE
effectsize::cohens_d(behavioral_prej ~ treatment, data = df)

# Visualize the results
group_mean_temp2 <- df %>% 
  group_by(treatment) %>% 
  summarize(xvalue = mean(behavioral_prej, na.rm = T))
p2 <- 
  ggplot(data = df, aes(x = behavioral_prej, group = treatment, fill = treatment)) +
  geom_density(alpha = .25) +
  geom_vline(data = group_mean_temp2,
             aes(xintercept = xvalue, color = treatment), 
             linetype = "dashed", linewidth = 0.6) +
  scale_color_manual(labels = c("Control", "Treatment"), values = c("#00008B", "#EFC000FF")) +
  scale_fill_manual(labels = c("Control", "Treatment"), values = c("#00008B", "#EFC000FF")) +
  xlab("Desire for Social Distance from Immigrants\n(0 = lowest; 100 = highest)") +
  ylab("Density") +
  ggtitle("Desire for Social Distance in Full Sample") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.95, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic"))
p2 <- p2 + 
  annotate(geom = "text", x = 75, y = 0.026, 
           label = paste("beta == ", round(beh_prej_ATE$coefficients[2], 2)),
           size = 4, family = "Times", parse = T) +
  annotate(geom = "text", x = 75, y = 0.0235, 
           label = paste("(SE == ", round(beh_prej_ATE$std.error[2], 2), ")"),
           size = 4, family = "Times", parse = T)
main_ATE <- plot_grid(p1, p2, labels = "AUTO", nrow = 1, label_fontfamily = "Times")
ggsave(file = "Figure 6.pdf", main_ATE, width = 9, height = 4.5)
# ggsave(file = "Figure 6.tif", main_ATE, dpi = 300, width = 9, height = 4.5)

# CATE by party ID
lm_robust(behavioral_prej ~ treatment + age + female + white + married +
            education + income + ideo,
          data = df,
          subset = (dem == 1))
effectsize::cohens_d(behavioral_prej ~ treatment, data = subset(df, dem == 1))

lm_robust(behavioral_prej ~ treatment + age + female + white + married +
            education + income + ideo,
          data = df,
          subset = (dem == 0))
effectsize::cohens_d(behavioral_prej ~ treatment, data = subset(df, dem == 0))

# HTE by seven-point party ID
out <- interflex(Y = "behavioral_prej", D = "treatment", X = "pid",
                 Z = c("age", "female", "white", "married", "education", "income", "ideo"), 
                 data = df, na.rm = TRUE,
                 estimator = "binning", vcov.type = "robust")
p2_HTE <- 
  plot(out, theme.bw = TRUE, cex.axis = 0.8, cex.lab = 0.8, bin.labs = TRUE,
       xlab = "Party Identification\n(0 = strong Democrat; 6 = strong Republican)", 
       ylab = "Marginal Effect on Desire for Social Distance",
       ylim = c(-15, 15)) +
  ggtitle("HTE on Desire for Social Distance by Party Identification") + 
  theme(plot.title = element_text(hjust = 0.85, size = 12, face = "italic", family = "Times"))
HTE_pid <- plot_grid(p1_HTE, p2_HTE, labels = "AUTO", nrow = 1, label_fontfamily = "Times")
ggsave(file = "Figure S16.pdf", HTE_pid, width = 9, height = 4.5)

# CATE by racial identity
lm_robust(behavioral_prej ~ treatment + age + female + married +
            education + income + ideo + pid,
          data = df,
          subset = (white == 1))
effectsize::cohens_d(behavioral_prej ~ treatment, data = subset(df, white == 1))

lm_robust(behavioral_prej ~ treatment + age + female + married +
            education + income + ideo + pid,
          data = df,
          subset = (white == 0))
effectsize::cohens_d(behavioral_prej ~ treatment, data = subset(df, white == 0))

### Figure S25 ----
## Support for immigration
# Estimate the ATE
ATE_immi_support <- lm_robust(immi_support ~ treatment + age + female + white + 
                                married + education + income + ideo + pid,
                              data = df)

# Estimate the effect size of the ATE
effectsize::cohens_d(immi_support ~ treatment, data = df)

# Visualize the ATE
summary_immi_support <- df %>% 
  group_by(treatment) %>% 
  do(tidy(lm_robust(immi_support ~ 1, data = .))) %>% 
  mutate(immi_support = estimate)
p <- 
  ggplot(summary_immi_support, aes(x = treatment, y = immi_support)) +
  geom_point(data = df, size = 1, alpha = .05, color = "grey",
             position = position_jitter(width = .1, height = .2, seed = 1)) +
  geom_point(size = 1.5) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), linewidth = .5, width = 0) +
  theme_classic() +
  scale_x_discrete(labels = c("0" = "Control", "1" = "Treatment")) +
  xlab("") +
  ylab("Preferred Number of Immigrants\n(0 = should decrease a lot; 4 = should increase a lot)") +
  coord_cartesian(ylim = c(-.25, 4.25)) +
  theme(text = element_text(family = "Times", size = 12),
        axis.text = element_text(family = "Times", size = 12))
p <- p + 
  annotate(geom = "text", x = 1.5, y = 3.1, 
           label = paste("beta == ", round(ATE_immi_support$coefficients[2], 2)),
           size = 4, family = "Times", parse = T) +
  annotate(geom = "text", x = 1.5, y = 2.9, 
           label = paste("(SE == ", round(ATE_immi_support$std.error[2], 2), ")"),
           size = 4, family = "Times", parse = T)
ggsave(file = "Figure S25.pdf", width = 5, height = 5)

# CATE by party ID
lm_robust(immi_support ~ treatment + age + female + white + married +
            education + income + ideo,
          data = df,
          subset = (dem == 1))
effectsize::cohens_d(immi_support ~ treatment, data = subset(df, dem == 1))

lm_robust(immi_support ~ treatment + age + female + white + married +
            education + income + ideo,
          data = df,
          subset = (dem == 0))
effectsize::cohens_d(immi_support ~ treatment, data = subset(df, dem == 0))

# HTE by seven-point party ID
out <- interflex(Y = "immi_support", D = "treatment", X = "pid",
                 Z = c("age", "female", "white", "married", "education", "income", "ideo"), 
                 data = df, na.rm = TRUE,
                 estimator = "binning", vcov.type = "robust")
p <- 
  plot(out, theme.bw = TRUE, cex.axis = 0.8, cex.lab = 0.8, bin.labs = TRUE,
       xlab = "Party Identification", ylab = "Marginal Effect on Immigration Support") +
  ggtitle("HTE by Party Identification") + 
  theme(plot.title = element_text(hjust = 0.5, size = 13, face = "italic", family = "Times"))

# CATE by racial identity
lm_robust(immi_support ~ treatment + age + female + married +
            education + income + ideo + pid,
          data = df,
          subset = (white == 1))
effectsize::cohens_d(immi_support ~ treatment, data = subset(df, white == 1))

lm_robust(immi_support ~ treatment + age + female + married +
            education + income + ideo + pid,
          data = df,
          subset = (white == 0))
effectsize::cohens_d(immi_support ~ treatment, data = subset(df, white == 0))

### Table S7 ----
## Identify respondents who believed that undocumented immigrants / legal immigrants had a higher crime rate than locals
df$misper_UI <- ifelse(df$pre_est_3 > df$pre_est_1, 1, 0)
df$misper_LI <- ifelse(df$pre_est_2 > df$pre_est_1, 1, 0)

## Analyze the correlations between misperceptions and affective prejudice / desire for social distance
mod1.1 <- lm_robust(affective_prej ~ treatment, 
                    data = df)
mod1.2 <- lm_robust(affective_prej ~ treatment +
                      age + female + white + married + education + income + ideo + pid, 
                    data = df)
mod1.3 <- lm_robust(affective_prej ~ treatment +
                      age + female + white + married + education + income + ideo + pid +
                      misper_LI + misper_UI, 
                    data = df)
mod1.4 <- lm_robust(affective_prej ~ treatment +
                      age + female + white + married + education + income + ideo + pid +
                      misper_LI + misper_UI, 
                    data = df, subset = attentive == 2)
mod2.1 <- lm_robust(behavioral_prej ~ treatment, 
                    data = df)
mod2.2 <- lm_robust(behavioral_prej ~ treatment + 
                      age + female + white + married + education + income + ideo + pid, 
                    data = df)
mod2.3 <- lm_robust(behavioral_prej ~ treatment + 
                      age + female + white + married + education + income + ideo + pid +
                      misper_LI + misper_UI, 
                    data = df)
mod2.4 <- lm_robust(behavioral_prej ~ treatment + 
                      age + female + white + married + education + income + ideo + pid +
                      misper_LI + misper_UI, 
                    data = df, subset = attentive == 2)

## Tabularize the results
texreg(list(mod1.1, mod1.2, mod1.3, mod1.4, mod2.1, mod2.2, mod2.3, mod2.4),
       custom.coef.names = c("Intercept", "Treatment (= 1)", 
                             "Age", "Female (= 1)", "White (= 1)", "Married (= 1)", 
                             "Education (5 = highest)", "Income (12 = highest)", 
                             "Self-Reported Ideology", "Party Identification",
                             "Misperception about LI (= 1)", "Misperception about UI (= 1)"),
       custom.header = list("Dependent Variable" = 1:8),
       custom.model.names = c("Affective", "Affective", "Affective", "Affective",
                              "Behavioral", "Behavioral", "Behavioral", "Behavioral"),
       stars = numeric(0),
       include.adjrs = TRUE, include.ci = FALSE, include.rmse = FALSE, 
       caption = "Misperceptions about Immigrant Crime, Factual Corrections, and Prejudice against Immigrants")

### Figure S17 ----
## Standardize the outcome variables
df <- df %>% mutate(
  affective_prej_z = scale(affective_prej),
  behavioral_prej_z = scale(behavioral_prej),
  immi_support_z = scale(immi_support)
)

## Create an empty data frame to store the estimates
df_race <- as.data.frame(matrix(NA, nrow = 6, ncol = 5))
df_race <- df_race %>% 
  rename(race = V1, dv = V2, cate = V3, lower_ci = V4, upper_ci = V5)
df_race$race <- c("Nonwhite", "White", "Nonwhite", "White", "Nonwhite", "White")
df_race$dv <- c(1, 1, 2, 2, 3, 3)

## DV = affective prejudice
temp1 <- lm_robust(affective_prej_z ~ treatment * white + age + female + white + married +
                     education + income + ideo + pid,
                   data = df)
CATE_temp1_nonwhite <- glht(temp1, linfct = c("treatment1 = 0"))
CATE_temp1_white <- glht(temp1, linfct = c("treatment1 + treatment1:white = 0"))
df_race[1, 3:5] <- confint(CATE_temp1_nonwhite)$confint
df_race[2, 3:5] <- confint(CATE_temp1_white)$confint

## DV = desire for social distance
temp2 <- lm_robust(behavioral_prej_z ~ treatment * white + age + female + white + married +
                     education + income + ideo + pid,
                   data = df)
CATE_temp2_nonwhite <- glht(temp2, linfct = c("treatment1 = 0"))
CATE_temp2_white <- glht(temp2, linfct = c("treatment1 + treatment1:white = 0"))
df_race[3, 3:5] <- confint(CATE_temp2_nonwhite)$confint
df_race[4, 3:5] <- confint(CATE_temp2_white)$confint

## DV = immigration support
temp3 <- lm_robust(immi_support_z ~ treatment * white + age + female + white + married +
                     education + income + ideo + pid,
                   data = df)
CATE_temp3_nonwhite <- glht(temp3, linfct = c("treatment1 = 0"))
CATE_temp3_white <- glht(temp3, linfct = c("treatment1 + treatment1:white = 0"))
df_race[5, 3:5] <- confint(CATE_temp3_nonwhite)$confint
df_race[6, 3:5] <- confint(CATE_temp3_white)$confint

## Visualize the CATEs by racial identity and by dependent variable (all results)
df_race$dv <- factor(df_race$dv,
                     levels = c(1:3),
                     labels = c("Affective Prejudice", "Desire for Social Distance", "Immigration Support"))
p <- 
  ggplot(data = df_race, aes(x = dv, y = cate, color = race, shape = race)) +
  geom_point(position = position_dodge(.25), size = 3) +
  scale_color_manual(values = c("grey25", "grey75")) +
  scale_shape_manual(values = c(16, 17)) +
  geom_errorbar(width = 0, aes(ymin = lower_ci, ymax = upper_ci), 
                position = position_dodge(.25)) +
  xlab("") + 
  ylab("Conditional Average Treatment Effect") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.99, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  coord_cartesian(ylim = c(-.5, .5))
ggsave(file = "Figure S17.pdf", width = 8, height = 4)

### Figure S8 ----
## % of respondents disagreeing that "immigrants' crime rate > locals' crime rate" in full sample
df %>% group_by(treatment) %>% summarize(n = n(), xvalue = mean(mani_check, na.rm = T))
lm_robust(mani_check ~ treatment + age + female + white + married +
            education + income + ideo + pid,
          data = df)

## % of respondents disagreeing that "immigrants' crime rate > locals' crime rate" by party ID
df_manip <- as.data.frame(matrix(NA, nrow = 4, ncol = 5))
df_manip <- df_manip %>% 
  rename(gop = V1, group = V2, percent = V3, lower_ci = V4, upper_ci = V5)
df_manip$gop <- c("Republican", "Republican", "Non-Republican", "Non-Republican")
df_manip$group <- c("Control", "Treatment", "Control", "Treatment")
temp1 <- prop.test(x = sum(df$mani_check[df$treatment == 0 & df$gop == 1], na.rm = T),
                   n = sum(!is.na(df$mani_check[df$treatment == 0 & df$gop == 1])),
                   correct = F)
temp2 <- prop.test(x = sum(df$mani_check[df$treatment == 1 & df$gop == 1], na.rm = T),
                   n = sum(!is.na(df$mani_check[df$treatment == 1 & df$gop == 1])),
                   correct = F)
temp3 <- prop.test(x = sum(df$mani_check[df$treatment == 0 & df$gop == 0], na.rm = T),
                   n = sum(!is.na(df$mani_check[df$treatment == 0 & df$gop == 0])),
                   correct = F)
temp4 <- prop.test(x = sum(df$mani_check[df$treatment == 1 & df$gop == 0], na.rm = T),
                   n = sum(!is.na(df$mani_check[df$treatment == 1 & df$gop == 0])),
                   correct = F)
df_manip[1, 3] <- temp1$estimate
df_manip[2, 3] <- temp2$estimate
df_manip[3, 3] <- temp3$estimate
df_manip[4, 3] <- temp4$estimate
df_manip[1, 4:5] <- temp1$conf.int
df_manip[2, 4:5] <- temp2$conf.int
df_manip[3, 4:5] <- temp3$conf.int
df_manip[4, 4:5] <- temp4$conf.int
df_manip$gop <- factor(df_manip$gop, levels = c("Republican", "Non-Republican"))
p_manip_pid <- 
  ggplot(data = df_manip, aes(x = gop, y = percent * 100, fill = group)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "black", alpha = .5) +
  scale_x_discrete(breaks = c("Republican", "Non-Republican")) +
  scale_fill_manual(values = c("grey15", "grey85")) +
  geom_errorbar(width = .2, aes(ymin = lower_ci * 100, ymax = upper_ci * 100), 
                position = position_dodge(.9)) +
  xlab("") + 
  ylab("Percent of Respondents Disagreeing that\nImmigrants' Crime Rate > Locals' Crime Rate") +
  ggtitle("Direct Manipulation Check by Partisanship") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.3, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic")) +
  coord_cartesian(ylim = c(0, 100))

## % of respondents disagreeing that "immigrants' crime rate > locals' crime rate" by racial identity
df_manip <- as.data.frame(matrix(NA, nrow = 4, ncol = 5))
df_manip <- df_manip %>% 
  rename(gop = V1, group = V2, percent = V3, lower_ci = V4, upper_ci = V5)
df_manip$gop <- c("White", "White", "Nonwhite", "Nonwhite")
df_manip$group <- c("Control", "Treatment", "Control", "Treatment")
temp1 <- prop.test(x = sum(df$mani_check[df$treatment == 0 & df$white == 1], na.rm = T),
                   n = sum(!is.na(df$mani_check[df$treatment == 0 & df$white == 1])),
                   correct = F)
temp2 <- prop.test(x = sum(df$mani_check[df$treatment == 1 & df$white == 1], na.rm = T),
                   n = sum(!is.na(df$mani_check[df$treatment == 1 & df$white == 1])),
                   correct = F)
temp3 <- prop.test(x = sum(df$mani_check[df$treatment == 0 & df$white == 0], na.rm = T),
                   n = sum(!is.na(df$mani_check[df$treatment == 0 & df$white == 0])),
                   correct = F)
temp4 <- prop.test(x = sum(df$mani_check[df$treatment == 1 & df$white == 0], na.rm = T),
                   n = sum(!is.na(df$mani_check[df$treatment == 1 & df$white == 0])),
                   correct = F)
df_manip[1, 3] <- temp1$estimate
df_manip[2, 3] <- temp2$estimate
df_manip[3, 3] <- temp3$estimate
df_manip[4, 3] <- temp4$estimate
df_manip[1, 4:5] <- temp1$conf.int
df_manip[2, 4:5] <- temp2$conf.int
df_manip[3, 4:5] <- temp3$conf.int
df_manip[4, 4:5] <- temp4$conf.int
df_manip$gop <- factor(df_manip$gop, levels = c("White", "Nonwhite"))
p_manip_race <- 
  ggplot(data = df_manip, aes(x = gop, y = percent * 100, fill = group)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "black", alpha = .5) +
  scale_x_discrete(breaks = c("White", "Nonwhite")) +
  scale_fill_manual(values = c("grey15", "grey85")) +
  geom_errorbar(width = .2, aes(ymin = lower_ci * 100, ymax = upper_ci * 100), 
                position = position_dodge(.9)) +
  xlab("") + 
  ylab("Percent of Respondents Disagreeing that\nImmigrants' Crime Rate > Locals' Crime Rate") +
  ggtitle("Direct Manipulation Check by Racial Identity") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.3, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic")) +
  coord_cartesian(ylim = c(0, 100))
manip_check <- plot_grid(p_manip_pid, p_manip_race, labels = "AUTO", nrow = 1, label_fontfamily = "Times")
ggsave(file = "Figure S8.pdf", manip_check, width = 9, height = 4.5)

### Figure S9 ----
## Perceived effect of illegal immigration on crime rates in full sample
df %>% group_by(treatment) %>% summarize(n = n(), xvalue = mean(perceived_crime, na.rm = T))
lm_robust(perceived_crime ~ treatment + age + female + white + married +
            education + income + ideo + pid,
          data = df)

## Perceived effect of illegal immigration on crime rates by party ID
summary_perceived_crime <- df %>% 
  group_by(treatment, gop) %>% 
  do(tidy(lm_robust(perceived_crime ~ 1, data = .))) %>% 
  mutate(perceived_crime = estimate) %>% 
  na.omit()
summary_perceived_crime$gop <- factor(summary_perceived_crime$gop, levels = c(1, 0), 
                                      labels = c("Republican", "Non-Republican"))
summary_perceived_crime$treatment <- factor(summary_perceived_crime$treatment, levels = c(0, 1), 
                                            labels = c("Control", "Treatment"))
p_manip_pid_2 <- 
  ggplot(data = summary_perceived_crime, aes(x = gop, y = perceived_crime, fill = treatment)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "black", alpha = .5) +
  scale_x_discrete(breaks = c("Republican", "Non-Republican")) +
  scale_fill_manual(values = c("grey15", "grey85")) +
  geom_errorbar(width = .2, aes(ymin = conf.low, ymax = conf.high), 
                position = position_dodge(.9)) +
  xlab("") + 
  ylab("Perceived Effect of Illegal Immigration on Crime\n(0 = decrease a lot; 6 = increase a lot)") +
  ggtitle("Posterior Beliefs by Partisanship") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.3, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic")) +
  coord_cartesian(ylim = c(0, 6))

## Perceived effect of illegal immigration on crime rates by racial identity
summary_perceived_crime <- df %>% 
  group_by(treatment, white) %>% 
  do(tidy(lm_robust(perceived_crime ~ 1, data = .))) %>% 
  mutate(perceived_crime = estimate) %>% 
  na.omit()
summary_perceived_crime$white <- factor(summary_perceived_crime$white, levels = c(1, 0), 
                                        labels = c("White", "Nonwhite"))
summary_perceived_crime$treatment <- factor(summary_perceived_crime$treatment, levels = c(0, 1), 
                                            labels = c("Control", "Treatment"))
p_manip_race_2 <- 
  ggplot(data = summary_perceived_crime, aes(x = white, y = perceived_crime, fill = treatment)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "black", alpha = .5) +
  scale_x_discrete(breaks = c("White", "Nonwhite")) +
  scale_fill_manual(values = c("grey15", "grey85")) +
  geom_errorbar(width = .2, aes(ymin = conf.low, ymax = conf.high), 
                position = position_dodge(.9)) +
  xlab("") + 
  ylab("Perceived Effect of Illegal Immigration on Crime\n(0 = Decrease a Lot; 6 = Increase a Lot)") +
  ggtitle("Posterior Beliefs by Racial Identity") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.3, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic")) +
  coord_cartesian(ylim = c(0, 6))
manip_check <- plot_grid(p_manip_pid_2, p_manip_race_2, labels = "AUTO", nrow = 1, label_fontfamily = "Times")
ggsave(file = "Figure S9.pdf", manip_check, width = 9, height = 4.5)

### Figure S26 ----
## Perceived credibility in full sample
df %>% group_by(treatment) %>% summarize(n = n(), xvalue = mean(cred, na.rm = T))
lm_robust(cred ~ treatment + age + female + white + married +
            education + income + ideo + pid,
          data = df)

# Republicans only
lm_robust(cred ~ treatment + age + female + white + married +
            education + income + ideo + pid,
          data = df, subset = gop == 1)

# Non-Republicans only
lm_robust(cred ~ treatment + age + female + white + married +
            education + income + ideo + pid,
          data = df, subset = gop == 0)

## Perceived credibility of the vignette by partisanship
summary_credibility <- df %>% 
  group_by(treatment, gop) %>% 
  do(tidy(lm_robust(cred ~ 1, data = .))) %>% 
  mutate(cred = estimate) %>% 
  na.omit()
summary_credibility$gop <- factor(summary_credibility$gop, levels = c(1, 0), 
                                    labels = c("Republican", "Non-Republican"))
summary_credibility$treatment <- factor(summary_credibility$treatment, levels = c(0, 1), 
                                        labels = c("Control", "Treatment"))
p_cred_pid <- 
  ggplot(data = summary_credibility, aes(x = gop, y = cred, fill = treatment)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "black", alpha = .5) +
  scale_x_discrete(breaks = c("Republican", "Non-Republican")) +
  scale_fill_manual(values = c("grey15", "grey85")) +
  geom_errorbar(width = .2, aes(ymin = conf.low, ymax = conf.high), 
                position = position_dodge(.9)) +
  xlab("") + 
  ylab("Perceived Credibility of Information\n(0 = very untrustworthy; 4 = very trustworthy)") +
  ggtitle("Perceived Credibility by Partisanship") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.3, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic")) +
  coord_cartesian(ylim = c(0, 4))

## Perceived credibility of the vignette by racial identity
summary_credibility <- df %>% 
  group_by(treatment, white) %>% 
  do(tidy(lm_robust(cred ~ 1, data = .))) %>% 
  mutate(cred = estimate) %>% 
  na.omit()
summary_credibility$white <- factor(summary_credibility$white, levels = c(1, 0), 
                                    labels = c("White", "Nonwhite"))
summary_credibility$treatment <- factor(summary_credibility$treatment, levels = c(0, 1), 
                                        labels = c("Control", "Treatment"))
p_cred_race <- 
  ggplot(data = summary_credibility, aes(x = white, y = cred, fill = treatment)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "black", alpha = .5) +
  scale_x_discrete(breaks = c("White", "Nonwhite")) +
  scale_fill_manual(values = c("grey15", "grey85")) +
  geom_errorbar(width = .2, aes(ymin = conf.low, ymax = conf.high), 
                position = position_dodge(.9)) +
  xlab("") + 
  ylab("Perceived Credibility of Information\n(0 = Very Untrustworthy; 4 = Very Trustworthy)") +
  ggtitle("Perceived Credibility by Racial Identity") +
  theme_classic() +
  theme(text = element_text(family = "Times", size = 13),
        axis.text = element_text(color = "black", size = 13),
        legend.justification = c(1, 1), legend.position = c(.3, .99),
        legend.box.background = element_rect(color = "black"), 
        legend.key.size = unit(1, "line"),
        legend.direction = "vertical",
        legend.margin = margin(t = 0.2, r = 0.2, b = 0.2, l = 0.2, unit = "cm"),
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5, size = 13, face = "italic")) +
  coord_cartesian(ylim = c(0, 4))
manip_check <- plot_grid(p_cred_pid, p_cred_race, labels = "AUTO", nrow = 1, label_fontfamily = "Times")
ggsave(file = "Figure S26.pdf", manip_check, width = 9, height = 4.5)

### Figure S18 ----
## HTE on affective prejudice against legal immigrants
out1 <- interflex(Y = "affective_prej", D = "treatment", X = "overest_LI",
                  Z = c("age", "female", "white", "married", "education", "income", "ideo", "pid"), 
                  data = df, na.rm = TRUE,
                  estimator = "kernel", vcov.type = "robust")
p1_HTE <- 
  plot(out1, theme.bw = TRUE, cex.axis = 0.8, cex.lab = 0.8, bin.labs = TRUE,
       xlab = "Level of Overestimation", 
       ylab = "Marginal Effect on Affective Prejudice",
       ylim = c(-10, 5)) +
  ggtitle("HTE on Affective Prejudice by Level of\nOverestimation against Legal Immigrants") + 
  theme(plot.title = element_text(hjust = 0.6, size = 13, face = "italic", family = "Times"))

## Desire for social distance from legal immigrants
out2 <- interflex(Y = "behavioral_prej", D = "treatment", X = "overest_LI",
                  Z = c("age", "female", "white", "married", "education", "income", "ideo", "pid"), 
                  data = df, na.rm = TRUE,
                  estimator = "kernel", vcov.type = "robust")
p2_HTE <- 
  plot(out2, theme.bw = TRUE, cex.axis = 0.8, cex.lab = 0.8, bin.labs = TRUE,
       xlab = "Level of Overestimation", 
       ylab = "Marginal Effect on Desire for Social Distance",
       ylim = c(-10, 5)) +
  ggtitle("HTE on Desire for Social Distance by Level of\nOverestimation against Legal Immigrants") + 
  theme(plot.title = element_text(hjust = 0.6, size = 13, face = "italic", family = "Times"))

## Affective prejudice against undocumented immigrants
out3 <- interflex(Y = "affective_prej", D = "treatment", X = "overest_UI",
                  Z = c("age", "female", "white", "married", "education", "income", "ideo", "pid"), 
                  data = df, na.rm = TRUE,
                  estimator = "kernel", vcov.type = "robust")
p3_HTE <- 
  plot(out3, theme.bw = TRUE, cex.axis = 0.8, cex.lab = 0.8, bin.labs = TRUE,
       xlab = "Level of Overestimation", 
       ylab = "Marginal Effect on Affective Prejudice",
       ylim = c(-10, 5)) +
  ggtitle("HTE on Affective Prejudice by Level of\nOverestimation against Undocumented Immigrants") + 
  theme(plot.title = element_text(hjust = 0.6, size = 13, face = "italic", family = "Times"))

## Desire for social distance from undocumented immigrants
out4 <- interflex(Y = "behavioral_prej", D = "treatment", X = "overest_UI",
                  Z = c("age", "female", "white", "married", "education", "income", "ideo", "pid"), 
                  data = df, na.rm = TRUE,
                  estimator = "kernel", vcov.type = "robust")
p4_HTE <- 
  plot(out4, theme.bw = TRUE, cex.axis = 0.8, cex.lab = 0.8, bin.labs = TRUE,
       xlab = "Level of Overestimation", 
       ylab = "Marginal Effect on Desire for Social Distance",
       ylim = c(-10, 5)) +
  ggtitle("HTE on Desire for Social Distance by Level of\nOverestimation against Undocumented Immigrants") + 
  theme(plot.title = element_text(hjust = 0.6, size = 13, face = "italic", family = "Times"))
HTE_overest <- plot_grid(p1_HTE, p2_HTE, p3_HTE, p4_HTE, labels = "AUTO", nrow = 2, label_fontfamily = "Times")
ggsave(file = "Figure S18.pdf", HTE_overest, width = 10, height = 9)

### Table S10 ----
## CACE on affective prejudice
cace_1 <- ivreg(affective_prej ~ mani_check | treatment, data = df)
cace_2 <- ivreg(affective_prej ~ mani_check + age + female + white + married +
                  education + income + ideo + pid | 
                  treatment + age + female + white + married +
                  education + income + ideo + pid, data = df)

## CACE on desire for social distance
cace_3 <- ivreg(behavioral_prej ~ mani_check | treatment, data = df)
cace_4 <- ivreg(behavioral_prej ~ mani_check + age + female + white + married +
                  education + income + ideo + pid | 
                  treatment + age + female + white + married +
                  education + income + ideo + pid, data = df)

## Summarize the results
m_list <- list(IV1 = cace_1, IV2 = cace_2, IV3 = cace_3, IV4 = cace_4)
msummary(m_list, output = "markdown")
msummary(m_list, output = "latex", fmt = fmt_decimal(digits = 2))

## Calculate the effect size in terms of Cohen's d
# Affective prejudice
iv_coef_1 <- summary(cace_1)$coefficients["mani_check", "Estimate"]
iv_se_1 <- summary(cace_1)$coefficients["mani_check", "Std. Error"]
sd_affective_prej <- sd(df$affective_prej, na.rm = TRUE)
cohen_d_1 <- iv_coef_1 / sd_affective_prej
cohen_d_se_1 <- iv_se_1 / sd_affective_prej
lower_1 <- cohen_d_1 - 1.96 * cohen_d_se_1
upper_1 <- cohen_d_1 + 1.96 * cohen_d_se_1
cat("Cohen's d:", cohen_d_1, "\n")
cat("95% CI: [", lower_1, ", ", upper_1, "]\n")

# Desire for social distance
iv_coef_3 <- summary(cace_3)$coefficients["mani_check", "Estimate"]
iv_se_3 <- summary(cace_3)$coefficients["mani_check", "Std. Error"]
sd_behavioral_prej <- sd(df$behavioral_prej, na.rm = TRUE)
cohen_d_3 <- iv_coef_3 / sd_behavioral_prej
cohen_d_se_3 <- iv_se_3 / sd_behavioral_prej
lower_3 <- cohen_d_3 - 1.96 * cohen_d_se_3
upper_3 <- cohen_d_3 + 1.96 * cohen_d_se_3
cat("Cohen's d:", cohen_d_3, "\n")
cat("95% CI: [", lower_3, ", ", upper_3, "]\n")

## Support for immigration
cace_5 <- ivreg(immi_support ~ mani_check | treatment, data = df)
msummary(cace_5)

iv_coef_5 <- summary(cace_5)$coefficients["mani_check", "Estimate"]
iv_se_5 <- summary(cace_5)$coefficients["mani_check", "Std. Error"]
sd_immi_support <- sd(df$immi_support, na.rm = TRUE)
cohen_d_5 <- iv_coef_5 / sd_immi_support
cohen_d_se_5 <- iv_se_5 / sd_immi_support
lower_5 <- cohen_d_5 - 1.96 * cohen_d_se_5
upper_5 <- cohen_d_5 + 1.96 * cohen_d_se_5
cat("Cohen's d:", cohen_d_5, "\n")
cat("95% CI: [", lower_5, ", ", upper_5, "]\n")
